package cn.zyl.demo.common.gif;

/** NeuQuant 神经网络量化算法<br/>
 * ------------------------------------------<br/>
 *
 * 版权所有（c）1994年安东尼·德克尔<br/>
 * NeuQuant 神经网络中的“网络：在计算神经系统”，由安东尼·德克尔量化算法，1994年见“获得最佳的色彩量化Kohonen神经网络”卷。 5（1994），第351-367。 该算法的讨论。<br/>
 * 任何一方从笔者获得这些文件的副本，直接或间接地，授予，免费的，完整的和不受限制的不可撤销的，全球范围内，缴纳的，免版税的，非独家的权利和许可，
 * 以经营本软件文档文件（“软件”），包括但不限于使用权，复制，修改，合并，发布，分发，授权和/或软件的销售拷贝，并允许谁收到副本来自任何一方的人员这样做，唯一的要求是本版权声明保持不变<br/>
 * Ported to Java 12/00 K Wei ner
 * @author wsz
 * @version v1.0
 * @date 2020-08-19 <a href="https://blog.csdn.net/u012726702/article/details/52015616">源码</a>
 */
public class NeuQuant {
    /** 使用的颜色数 */
    protected static final int NET_SIZE = 256;
    /** 接近500的四个素数-假设没有图像的长度如此之大,它可以被所有四个素数整除 */
    protected static final int PRIME_1 = 499;
    /** 接近500的素数 */
    protected static final int PRIME_2 = 491;
    /** 接近500的素数 */
    protected static final int PRIME_3 = 487;
    /** 接近500的素数 */
    protected static final int PRIME_4 = 503;
    /** 输入图像的最小尺寸*/
    protected static final int MIN_PICTURE_BYTES = (3 * PRIME_4);
    /* ------程序骨架 ---------------- */
    /** 网络定义 */
    protected static final int MAX_NET_POS = (NET_SIZE - 1);
    /** 颜色值的偏差 */
    protected static final int NET_BIAS_SHIFT = 4;
    /** 没有学习周期 */
    protected static final int NO_CYCLES = 100;
    /** 定义频率和偏差 */
    protected static final int INT_BIAS_SHIFT = 16;
    /** 分数偏差 */
    protected static final int INT_BIAS = (1 << INT_BIAS_SHIFT);
    protected static final int GAMMA_SHIFT = 10;
    /** gamma = 1024 */
    protected static final int GAMMA = (1 << GAMMA_SHIFT);
    protected static final int BETA_SHIFT = 10;
    /** beta = 1/1024 */
    protected static final int BETA = (INT_BIAS >> BETA_SHIFT);
    protected static final int BETA_GAMMA = (INT_BIAS << (GAMMA_SHIFT - BETA_SHIFT));
    /** def s 用于减小半径系数 */
    protected static final int INIT_RAD = (NET_SIZE >> 3);
    /** 对于256列cols，半径开始 */
    protected static final int RADIUS_BIAS_SHIFT = 6;
    /** 在32.0偏移6位 */
    protected static final int RADIUS_BIAS = (1 << RADIUS_BIAS_SHIFT);
    /** 并减少一个 */
    protected static final int INIT_RADIUS = (INIT_RAD * RADIUS_BIAS);
    /** 每个周期的1/30 */
    protected static final int RADIUS_DEC = 30;
    /** def s 用于降低alpha因子 */
    protected static final int ALPHA_BIAS_SHIFT = 10;
    /** alpha从1.0开始 */
    protected static final int INIT_ALPHA = (1 << ALPHA_BIAS_SHIFT);
    /** 偏差10位 */
    protected int alphaDec;
    /** rad偏置和alpha rad偏置用于rad功率计算 */
    protected static final int RAD_BIAS_SHIFT = 8;
    protected static final int RAD_BIAS = (1 << RAD_BIAS_SHIFT);
    protected static final int ALPHA_RAD_B_SHIFT = (ALPHA_BIAS_SHIFT + RAD_BIAS_SHIFT);
    protected static final int ALPHA_RAD_BIAS = (1 << ALPHA_RAD_B_SHIFT);

    /* ------类型和全局变量------ */
    /** 输入图像本身 */
    protected byte[] thePicture;
    /** 长度计数 = H*W*3 */
    protected int lengthCount;
    /** 采样系数 1..30 */
    protected int sampleFac;
    /** 类型def int 像素[4]; BGRc */
    protected int[][] network;
    /** 网络本身 - [netSize][4] */
    protected int[] netIndex = new int[256];
    /** 用于网络查找 - really 256 */
    protected int[] bias = new int[NET_SIZE];
    /** 偏差和频率数组用于学习 */
    protected int[] freq = new int[NET_SIZE];
    /** 预先计算的功率 */
    protected int[] radPower = new int[INIT_RAD];
    /** 初始化（0,0,0）到（255,255,255）范围内的网络并设置参数
     *<pre>
     *     [在1..30范围内选择样本fac]
     *     [从输入文件中读取图像]
     *     pic = (unsigned char*) malloc(3*width*height);
     *     initNet(pic,3*width*height,sampleFac);
     *     learn();
     *     unBiasNet();
     *     [使用写输出图像标题writeColourMap(f)]
     *     inxBuild();
     *     使用写输出图像 inxSearch(b,g,r)
     *</pre>
     */
    public NeuQuant(byte[] thePic, int len, int sample) {
        int i;
        int[] p;
        thePicture = thePic;
        lengthCount = len;
        sampleFac = sample;
        network = new int[NET_SIZE][];
        for (i = 0; i < NET_SIZE; i++) {
            network[i] = new int[4];
            p = network[i];
            p[0] = p[1] = p[2] = (i << (NET_BIAS_SHIFT + 8)) / NET_SIZE;
            /* 1/netSize */
            freq[i] = INT_BIAS / NET_SIZE;
            bias[i] = 0;
        }
    }
    public byte[] colorMap() {
        byte[] map = new byte[3 * NET_SIZE];
        int[] index = new int[NET_SIZE];
        for (int i = 0; i < NET_SIZE; i++) {
            index[network[i][3]] = i;
        }
        int k = 0;
        for (int i = 0; i < NET_SIZE; i++) {
            int j = index[i];
            map[k++] = (byte) (network[j][0]);
            map[k++] = (byte) (network[j][1]);
            map[k++] = (byte) (network[j][2]);
        }
        return map;
    }
    /** 网络的插入类型和netIndex [0..255]的构建（取消偏置后要做）*/
    public void inxBuild() {
        int i, j, smallpos, smallval, previouscol = 0, startpos= 0;
        int[] p,q;
        for (i = 0; i < NET_SIZE; i++) {
            p = network[i];
            smallpos = i;
            /* g指数 */
            smallval = p[1];
            /* 在i..netSize-1中找到最小的 */
            for (j = i + 1; j < NET_SIZE; j++) {
                q = network[j];
                /* g指数*/
                if (q[1] < smallval) {
                    smallpos = j;
                    smallval = q[1];
                }
            }
            q = network[smallpos];
            /* 交换p（i）和q（smallPos）条目 */
            if (i != smallpos) {
                j = q[0];
                q[0] = p[0];
                p[0] = j;
                j = q[1];
                q[1] = p[1];
                p[1] = j;
                j = q[2];
                q[2] = p[2];
                p[2] = j;
                j = q[3];
                q[3] = p[3];
                p[3] = j;
            }
            /* smallVal条目现在位于位置i */
            if (smallval != previouscol) {
                netIndex[previouscol] = (startpos + i) >> 1;
                for (j = previouscol + 1; j < smallval; j++) {
                    netIndex[j] = i;
                }
                previouscol = smallval;
                startpos = i;
            }
        }
        netIndex[previouscol] = (startpos + MAX_NET_POS) >> 1;
        for (j = previouscol + 1; j < 256; j++) {
            /* really 256 */
            netIndex[j] = MAX_NET_POS;
        }
    }
    /** 主要学习循环 */
    public void learn() {
        int i, j, b, g, r;
        int radius, rad, alpha, step, delta, samPlePixels;
        byte[] p;
        int pix, lim;
        if (lengthCount < MIN_PICTURE_BYTES) {
            sampleFac = 1;
        }
        alphaDec = 30 + ((sampleFac - 1) / 3);
        p = thePicture;
        pix = 0;
        lim = lengthCount;
        //样本像素
        samPlePixels = lengthCount / (3 * sampleFac);
        delta = samPlePixels / NO_CYCLES;
        alpha = INIT_ALPHA;
        radius = INIT_RADIUS;
        rad = radius >> RADIUS_BIAS_SHIFT;
        if (rad <= 1) {
            rad = 0;
        }
        for (i = 0; i < rad; i++) {
            radPower[i] = alpha * (((rad * rad - i * i) * RAD_BIAS) / (rad * rad));
        }
        //printf(stderr,"beginning 1D learning: initial radius=%d/n", rad);
        if (lengthCount < MIN_PICTURE_BYTES) {
            step = 3;
        } else if ((lengthCount % PRIME_1) != 0) {
            step = 3 * PRIME_1;
        }else {
            if ((lengthCount % PRIME_2) != 0) {
                step = 3 * PRIME_2;
            }else {
                if ((lengthCount % PRIME_3) != 0) {
                    step = 3 * PRIME_3;
                }else {
                    step = 3 * PRIME_4;
                }
            }
        }
        i = 0;
        while (i < samPlePixels) {
            b = (p[pix] & 0xff) << NET_BIAS_SHIFT;
            g = (p[pix + 1] & 0xff) << NET_BIAS_SHIFT;
            r = (p[pix + 2] & 0xff) << NET_BIAS_SHIFT;
            j = contest(b, g, r);
            alterSingle(alpha, j, b, g, r);
            if (rad != 0) {
                /* 改变邻居 */
                alterNeigh(rad, j, b, g, r);
            }
            pix += step;
            if (pix >= lim) {
                pix -= lengthCount;
            }
            i++;
            if (delta == 0) {
                delta = 1;
            }
            if (i % delta == 0) {
                alpha -= alpha / alphaDec;
                radius -= radius / RADIUS_DEC;
                rad = radius >> RADIUS_BIAS_SHIFT;
                if (rad <= 1) {
                    rad = 0;
                }
                for (j = 0; j < rad; j++) {
                    radPower[j] =
                            alpha * (((rad * rad - j * j) * RAD_BIAS) / (rad * rad));
                }
            }
        }
        // printf(stderr,"finished 1D learning: final alpha=%f !/n",((float)alpha)/initAlpha);
    }
    /** 搜索BGR值0..255（网络无偏后）并返回颜色索引 */
    public int map(int b, int g, int r) {
        int[] p;
        /* g指数 ；从netIndex [g]开始并向外工作 ; dist ; a ; 最大可能的dist是256 * 3;best */
        int i= netIndex[g], j = i-1, dist, a, bestd = 1000,best = -1;

        while ((i < NET_SIZE) || (j >= 0)) {
            if (i < NET_SIZE) {
                p = network[i];
                /* inx key */
                dist = p[1] - g;
                if (dist >= bestd) {
                    /* 停止迭代 */
                    i = NET_SIZE;
                }else {
                    i++;
                    if (dist < 0) {
                        dist = -dist;
                    }
                    a = p[0] - b;
                    if (a < 0) {
                        a = -a;
                    }
                    dist += a;
                    if (dist < bestd) {
                        a = p[2] - r;
                        if (a < 0) {
                            a = -a;
                        }
                        dist += a;
                        if (dist < bestd) {
                            bestd = dist;
                            best = p[3];
                        }
                    }
                }
            }
            if (j >= 0) {
                p = network[j];
                /* inx key - reverse dif */
                dist = g - p[1];
                if (dist >= bestd) {
                    /* 停止迭代 */
                    j = -1;
                }else {
                    j--;
                    if (dist < 0) {
                        dist = -dist;
                    }
                    a = p[0] - b;
                    if (a < 0) {
                        a = -a;
                    }
                    dist += a;
                    if (dist < bestd) {
                        a = p[2] - r;
                        if (a < 0) {
                            a = -a;
                        }
                        dist += a;
                        if (dist < bestd) {
                            bestd = dist;
                            best = p[3];
                        }
                    }
                }
            }
        }
        return (best);
    }
    public byte[] process() {
        learn();
        unBiasNet();
        inxBuild();
        return colorMap();
    }
    /** 无偏网络给出字节值0..255并记录位置i以准备进行排序 */
    public void unBiasNet() {
        for (int i = 0; i < NET_SIZE; i++) {
            network[i][0] >>= NET_BIAS_SHIFT;
            network[i][1] >>= NET_BIAS_SHIFT;
            network[i][2] >>= NET_BIAS_SHIFT;
            /* 记录颜色 */
            network[i][3] = i;
        }
    }
    /** 通过radPower [| i-j |]中预先计算的alpha *（1-（（i-j）^ 2 / [r] ^ 2））移动相邻的神经元 */
    protected void alterNeigh(int rad, int i, int b, int g, int r) {
        int j = i + 1, k = i - 1, lo = i - rad, hi = i + rad, a, m = 1;
        int[] p;
        if (lo < -1) {
            lo = -1;
        }
        if (hi > NET_SIZE) {
            hi = NET_SIZE;
        }
        while ((j < hi) || (k > lo)) {
            a = radPower[m++];
            if (j < hi) {
                p = network[j++];
                try {
                    p[0] -= (a * (p[0] - b)) / ALPHA_RAD_BIAS;
                    p[1] -= (a * (p[1] - g)) / ALPHA_RAD_BIAS;
                    p[2] -= (a * (p[2] - r)) / ALPHA_RAD_BIAS;
                    // 防止1.3错误编译
                } catch (Exception e) {
                    e.printStackTrace();
                }
            }
            if (k > lo) {
                p = network[k--];
                try {
                    p[0] -= (a * (p[0] - b)) / ALPHA_RAD_BIAS;
                    p[1] -= (a * (p[1] - g)) / ALPHA_RAD_BIAS;
                    p[2] -= (a * (p[2] - r)) / ALPHA_RAD_BIAS;
                } catch (Exception e) {
                    e.printStackTrace();
                }
            }
        }
    }
    /** 将神经元i移至因因子α偏向（b，g，r） */
    protected void alterSingle(int alpha, int i, int b, int g, int r) {
        /* 改变命中神经元 */
        int[] n = network[i];
        n[0] -= (alpha * (n[0] - b)) / INIT_ALPHA;
        n[1] -= (alpha * (n[1] - g)) / INIT_ALPHA;
        n[2] -= (alpha * (n[2] - r)) / INIT_ALPHA;
    }
    /** 搜索偏差的BGR值 */
    protected int contest(int b, int g, int r) {
        /* 查找最近的神经元（最小距离）并更新频率 */
        /* 找到最佳神经元（最小距离偏差）并返回位置 */
        /* 对于经常选择的神经元，freq [i]高，bias [i]负 */
        /* bias[i] = gamma*((1/netsize)-freq[i]) */
        int i, dist, a, biasdist, betafreq;
        int bestpos, bestbiaspos, bestd, bestbiasd;
        int[] n;
        bestd = ~(1 << 31);
        bestbiasd = bestd;
        bestpos = -1;
        bestbiaspos = bestpos;
        for (i = 0; i < NET_SIZE; i++) {
            n = network[i];
            dist = n[0] - b;
            if (dist < 0) {
                dist = -dist;
            }
            a = n[1] - g;
            if (a < 0) {
                a = -a;
            }
            dist += a;
            a = n[2] - r;
            if (a < 0) {
                a = -a;
            }
            dist += a;
            if (dist < bestd) {
                bestd = dist;
                bestpos = i;
            }
            biasdist = dist - ((bias[i]) >> (INT_BIAS_SHIFT - NET_BIAS_SHIFT));
            if (biasdist < bestbiasd) {
                bestbiasd = biasdist;
                bestbiaspos = i;
            }
            betafreq = (freq[i] >> BETA_SHIFT);
            freq[i] -= betafreq;
            bias[i] += (betafreq << GAMMA_SHIFT);
        }
        freq[bestpos] += BETA;
        bias[bestpos] -= BETA_GAMMA;
        return (bestbiaspos);
    }
}
